#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <algorithm>
#include <TCanvas.h>
#include <TGraphErrors.h>
#include <TMultiGraph.h>
#include <TStyle.h>
#include <TH1.h>
#include <TFile.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TF1.h>
#include <TTree.h>
#include <TChain.h>
#include <TFileCollection.h>
#include <TGraph2D.h>
#include <TLatex.h>
#include <TMarker.h>
#include <math.h>
#include <TMath.h>
#include <TFileInfo.h>

///Igor's stuff
#include "TROOT.h"
#include "MVATrainer.h"

#ifdef MVA
#include "MVAComputer2.h"
#endif

#include <TLorentzVector.h>
#include "MathRoot.h"
#include "TInterpreter.h"
#include "MSSM_tools/XSCalculator.C"
#include "MistagProbFunctions.h"
#include <sstream>
#include <iostream>
//#include "InputVariablesHistograms.h"
//#define TESTIGORCODE


//#include "../interface/tiltt.h"

//#define SUSY
//#define MC
//#define DATA

#define COUT(a) cout<<"\n"<<gInterpreter->GetCurrentMacroName()<<" :: "<<a<<"\n"<<"\n"

#define SMOOTH


#undef PTMAX300

#ifdef MC
#define LUMIWEIGHTING
//#define HLTWEIGHTING
//#define bPARTONMATCH
#endif

#ifdef DATA
//#define MATCHING
#endif

#include "readhbb.h"

TCanvas* canvas;

const int nflav=3;
const int novl = 4;

 int nInputEvents = 0;
 int nentries=0;
class bTagEff {
public:
  std::string file;
  std::string onloff;
//   const int& nflav;
//   const int& novl;
  const std::string* flavlabel;
  const std::string* sovl;

  TFile* hfile;
  TH2D* hrel2D[nflav][novl];
  //TH2D*** hrel2D;
bool _doPtEta;


  bTagEff(){};
  bTagEff(const std::string& theFile,const std::string& theonloff,const std::string* theFlavlabel,
          const std::string* theOvlLabel ){


 _doPtEta=false;
#if !defined(MISTAGPROBFUN) 
_doPtEta=true;
#endif


 if (theonloff == "online") _doPtEta=true;

if (_doPtEta) {

    if (theonloff != "online" && theonloff != "offline") {
      std::cout << "bTagEff constructor: Bad online/offline mode " << theonloff << std::endl;
      return;
    }


    hfile = TFile::Open( theFile.c_str() );
    if (hfile == NULL) {
      std::cout << "bTagEff constructor: Cannot open file " << theFile << std::endl;
      return;
    }

    onloff = theonloff;
    flavlabel = theFlavlabel;
    sovl = theOvlLabel;

    for (int iflav=0; iflav<nflav; ++iflav) {
      for (int iovl=0; iovl<novl; ++iovl) {
        if (theonloff == "online") {
#ifdef SMOOTH
          string shname( Form("pteta_%s_%s_releff_div-smp",sovl[iovl].c_str(),flavlabel[iflav].c_str()) );
          //string shname( Form("pteta_%s_%s_releff_sm0",sovl[iovl].c_str(),flavlabel[iflav].c_str()) );
          std::cout << "smoothed eff from " << shname << std::endl;
#else
          string shname( Form("pteta_%s_%s_releff_div",sovl[iovl].c_str(),flavlabel[iflav].c_str()) );
#endif

          hrel2D[iflav][iovl] = (TH2D*) hfile->Get( shname.c_str() );
          if (hrel2D[iflav][iovl] == NULL) std::cout << "Histogram not found: "
                            << shname.c_str() << std::endl;
          std::cout << "Using efficiency histogram " << hrel2D[iflav][iovl]->GetName() << " nbinsx= " << hrel2D[iflav][iovl]->GetXaxis()->GetNbins()
                    << " nbinsy= " << hrel2D[iflav][iovl]->GetYaxis()->GetNbins() << std::endl;
        } else {
          hrel2D[iflav][iovl] = (TH2D*) hfile->Get( Form("pteta_%s_%s_eff_div",sovl[iovl].c_str(),flavlabel[iflav].c_str()) );
          if (hrel2D[iflav][iovl] == NULL) std::cout << "Histogram not found: "
                                                     << Form("pteta_%s_%s_eff_div",sovl[iovl].c_str(),flavlabel[iflav].c_str()) << std::endl;
        }
      }
    }

} else 	COUT("Btag is only pT dependent");

  }
  double eff(int iflav,const std::string& theOvlLabel,double pt,double eta) {


if (_doPtEta) {


    if (iflav<0 || iflav>=nflav) {
      std::cout << "eff: bad flavor " << iflav << std::endl;
      return -999999;
    }
    int iovl = -1;
    for (int jovl=0; jovl<novl; ++jovl) {
        if (sovl[jovl] == theOvlLabel) {
          iovl = jovl;
          break;
        }
    }

    if (iovl == -1) {
      std::cout << "Bad btag type requested " << theOvlLabel << std::endl;
      return -999999;
    }
#ifndef PTMAX300
    // to avoid sampling outside of histogram borders
    if (pt>=300) pt=299.9;
#endif

    // just to test
    if ( hrel2D[iflav][iovl]->Interpolate(pt,fabs(eta)) == 0) {
      std::cout << "eff() " << theOvlLabel
                << "  " << onloff
                << " iflav=" << iflav
                << " zero for " << pt << " " << eta << std::endl;
      return 0;
    }

    return hrel2D[iflav][iovl]->Interpolate(pt,fabs(eta));

} else {

  switch (iflav) {
case 0:
	return pMistagUds(pt);
        break;
case 1:
	return pMistagC(pt);
        break;
case 2:
	return pMistagB(pt);
        break;
}

	return -999999;
	}

  }
};






class TrigHistArray {
public:
  std::vector<TH1F*> histos;
  TH1F* histAllTrig;
  TH1F* histAllTrigWeighted;
  std::vector<std::string>* gtl;  // genericTriggerList
  unsigned int filterPattern;
  
  TrigHistArray() {
    std::cout << "Standard constructor called, this should not happen" << std::endl;
  }

  TrigHistArray(std::vector<std::string>* theGtl,std::vector<std::string>* theFilter,const char* genName,const char* genTitle,int nx,float xmin,float xmax) {
    gtl = theGtl;
    char ctn[1000];
    char ctt[1000];
    // first set the filter pattern
    filterPattern = 0;
    if (theFilter != NULL) {
      for (std::vector<std::string>::iterator fit=theFilter->begin(); fit != theFilter->end(); ++fit) {
	std::vector<std::string>::iterator tSlotBtag = std::find(theGtl->begin(), theGtl->end(), *fit);
	if (tSlotBtag != theGtl->end()) {
	  filterPattern = filterPattern | (1<<(tSlotBtag - theGtl->begin()));
	} else {
	  std::cout << "Filter trigger " << *fit << " not found in any slot" << std::endl;
	}
      }
    } else {
      filterPattern = ~0;  // set all filter bits to one
    }

    for (unsigned int ib=0; ib<gtl->size(); ++ib) {
      TH1F* theHist = 0;
      if (filterPattern & (1<<ib)) {
	sprintf(ctn,"%sTrig%d",genName,ib);
	sprintf(ctt,"%s Trigger %s",genTitle,(*gtl)[ib].c_str());
	theHist = new TH1F(ctn,ctt,nx,xmin,xmax);
      }
      histos.push_back(theHist);
    }
    // inclusive histogram
    sprintf(ctn,"%sAllTrig",genName);
    sprintf(ctt,"%s all triggers",genTitle);
    histAllTrig = new TH1F(ctn,ctt,nx,xmin,xmax);
    sprintf(ctn,"%sAllTrigWeighted",genName);
    sprintf(ctt,"%s all triggers weighted",genTitle);
    histAllTrigWeighted = new TH1F(ctn,ctt,nx,xmin,xmax);
  }
			     
  void fill(int theTrgAccept,float x,float weight=1) {
    for (unsigned int ib=0; ib<gtl->size(); ++ib) {
      if (theTrgAccept & filterPattern & (1<<ib)) {
	if (histos[ib] != NULL) {
	  histos[ib]->Fill(x,weight);
	} else {
	  std::cout << "TrigHistArray::Fill: bad histogram pointer " << std::endl;
	}
      }
    }
    histAllTrig->Fill(x);
    histAllTrigWeighted->Fill(x,weight);
  }
};


void getHbbMetadata(TFileCollection& fc,std::vector<std::string>& genericTriggerList,const float targetLumi,
		    float& lumiScaleFac,const bool doMC) {

  // extracts the following metadata from the Hbb analysis ntuples:
  //     list of generic triggers (= dictionary for bit number convention)
  //     total number of input events before any filtering
  //     generator cross sections and filter efficiency
  // computes lumi scaling factor corresponding to given target integrated luminosity

  lumiScaleFac = 1;
//  int nInputEvents = 0;
  // first get the generic trigger list from the first file
  if (fc.GetNFiles() <=0) {
    std::cout << "getHbbMetadata: file list empty" << std::endl;
    return;
  }

  bool gtlDone = false;
  TIter iter( (TCollection*) fc.GetList() );
  while (TFileInfo* tfi = (TFileInfo*) &(*(iter.Next()))) {
    std::cout << "Open file " << tfi->GetCurrentUrl()->GetFile() << std::endl;
    std::cout << "Open Url " << tfi->GetCurrentUrl()->GetUrl() << std::endl;
    TFile* fg = TFile::Open( tfi->GetCurrentUrl()->GetUrl() );

    if (doMC) {
      TH1F* eventCount = (TH1F*) fg->Get("InputEvents/EventCount");
      if (eventCount != NULL) {
	std::cout << "Adding to event count " << eventCount->GetEntries() << std::endl;
	nInputEvents += eventCount->GetEntries();
      } else {
	std::cout << "eventCount histogram not found" << std::endl;
      }
    }

    if (! gtlDone) {
      TH1F* gtlHist = (TH1F*) fg->Get("hbbanalysis/gtlHist");
      if (gtlHist != NULL) {
	for (int ibin=1; ibin<= gtlHist->GetXaxis()->GetNbins(); ++ibin) {
	  genericTriggerList.push_back(gtlHist->GetXaxis()->GetBinLabel( ibin ));
	  std::cout << "getHbbMetadata: add generic trigger " << gtlHist->GetXaxis()->GetBinLabel( ibin )
		    << std::endl;
	}
	gtlDone = true;
      } else {
	std::cout << "gtlHist histogram not found" << std::endl;
      }
    }
    if (! doMC) break;  // for real data, only need first file to extract GTL
  }

  if (doMC) std::cout << "Total count of input events: " << nInputEvents << std::endl;
  std::cout << "Number of generic triggers read: " << genericTriggerList.size() << std::endl;

  if (doMC) {
    // now chain the GenInfo tree to access cross section & filter efficiency
    TChain gInfo("hbbanalysis/GenInfo");
    gInfo.AddFileInfoList((TCollection*) fc.GetList());
    TTree* gInfoTree = &gInfo;
    int nGInfo = gInfoTree->GetEntries();
    std::cout << "Number of events in gInfoTree: " << nGInfo << std::endl;

    // set GenInfo branches manually
    gInfoTree->SetBranchAddress("internalXsec",&_internalXsec);
    gInfoTree->SetBranchAddress("externalXsecLO",&_externalXsecLO);
    gInfoTree->SetBranchAddress("externalXsecNLO",&_externalXsecNLO);
    gInfoTree->SetBranchAddress("filterEfficiency",&_filterEfficiency);

    double SInternalXsec = 0;
    double SExternalXsecLO = 0;
    double SExternalXsecNLO = 0;
    double SFilterEfficiency = 0;
    for (int iGInfo=0; iGInfo<nGInfo; ++iGInfo) {
      gInfoTree->GetEntry(iGInfo);
      SInternalXsec += _internalXsec;
      SExternalXsecLO += _externalXsecLO;
      SExternalXsecNLO += _externalXsecNLO;
      SFilterEfficiency += _filterEfficiency;
    }
  
    double InternalXsec = SInternalXsec / nGInfo;
//     double ExternalXsecLO = SExternalXsecLO / nGInfo;
//     double ExternalXsecNLO = SExternalXsecNLO / nGInfo;
    double FilterEfficiency = SFilterEfficiency / nGInfo;

    std::cout << "Generator information:" << std::endl
	      << "      InternalXsec      = " << SInternalXsec / nGInfo << " pb " << std::endl
	      << "      ExternalXsecLO    = " << SExternalXsecLO / nGInfo << " pb " << std::endl
	      << "      ExternalXsecNLO   = " << SExternalXsecNLO / nGInfo << " pb " << std::endl
	      << "      FilterEfficiency  = " << SFilterEfficiency / nGInfo << std::endl;

    // in case not filter has been applied
    if (FilterEfficiency <0) {
      std::cout << "GenInfo FilterEfficieny <0 ==> set to 1" << std::endl;
      FilterEfficiency = 1;
    }

    // compute the lumi scale factor
    double XSec = InternalXsec;

///Igor's code
#ifdef SUSY
#ifdef SUSYMASS
	Double_t _mass = SUSYMASS;
//	   XSec=XSCalculator(_mass,20.,"MSSM_tools/");
	   XSec=XSCalculator(120.,20.,"MSSM_tools/"); ///always use 120 GeV mass. So called democratic mixture
#else
	cout<<"That's SUSY sample with standard mass 120"<<endl;
	XSec=111e0;
#endif

#endif

    double integLumi = nInputEvents / (XSec * FilterEfficiency);
    lumiScaleFac = targetLumi / integLumi;

    std::cout << "Lumi scale factor computation based on internal cross section:" << std::endl
	      << "      Xsec              = " << XSec << " pb " << std::endl
	      << "      IntLumiMC         = " << integLumi << " pb-1 " << std::endl
	      << "      TargetLumi        = " << targetLumi << " pb-1 " << std::endl
	      << "      Lumi Scale Factor = " << lumiScaleFac << std::endl;
  }
}

int jetFlavorCode(const int iJet) {
  int theFlavCode = -1;
  int theFlavor = partonFlavorJet[iJet];
  if (hflContentJet[iJet] != 0) theFlavor = hflContentJet[iJet];
  switch (abs(theFlavor)) {
  case 0:
  case 1:
  case 2:
  case 3:
    theFlavCode = 0;
    break;
  case 4:
    theFlavCode = 1;
    break;
  case 5:
    theFlavCode = 2;
    break;
  case 21:
    theFlavCode = 0;
    break;
  default:
    std::cout << "bad flavor " << theFlavor << std::endl;
  }
  return theFlavCode;
}

///Igor's code
TString PrintFlavor(int theFlavCode)
{
TString res=" not-defined jet ";
  switch (theFlavCode) {
case 0:	
	res = " udsg jet ";
	break;
case 1:
	 res = " c jet ";
	break;
case 2:
	res = " b jet ";
	break;
}

return res;
}


int diJetFlavorCode(const int iJet0, const int iJet1) {
  // determine dijet flavor from MC info
  // we set a flavor pair code that does not distinguish ordering
  int theFlav[2];
  theFlav[0] = jetFlavorCode( iJet0 );
  theFlav[1] = jetFlavorCode( iJet1 );
  // now set the dijet flavor
  int theFcDijet = -1;
  int flavPatternDijet = 10*theFlav[0]+theFlav[1];
  switch (flavPatternDijet) {
  case 0:
    theFcDijet = 5;
    break;
  case 11:
    theFcDijet = 3;
    break;
  case 22:
    theFcDijet = 0;
    break;
  case 1:
  case 10:
    theFcDijet = 4;
    break;
  case 2:
  case 20:
    theFcDijet = 2;
    break;
  case 12:
  case 21:
    theFcDijet = 1;
    break;
  default:
    std::cout << "diJetFlavorCode: Bad flavor code " << theFlav[0] << " " << theFlav[1] << std::endl;
  }
  return theFcDijet;
}


void readtribn() {

  // here we require three btagged jets



///Igor's code
Int_t _nevt_before_mva_trig=0;
Int_t _nevt_after_mva_trig=0;


  canvas = new TCanvas ("cg1","mycanvas",10,10,800,600);
  // open an ntuple file
  std::cout << " starting..." << std::endl;

  std::vector<std::string> genericTriggerList;
  std::vector<std::string> triggerFilterList;

#ifdef MC
  bool _doMC = true;
#else
  bool _doMC = false;
#endif

//cout<<"doMC = "<<_doMC<<endl;


  // create the trigger filter list. Only histograms for these triggers will be created
#ifdef MC
#ifdef BENRICHED
  triggerFilterList.push_back("HLT_CentralJet46_BTagIP3D_CentralJet38_BTagIP3D");
  triggerFilterList.push_back("HLT_CentralJet46_CentralJet38_DiBTagIP3D_v1");
#else
  triggerFilterList.push_back("HLT_bbPhi_CentralJet46_CentralJet38_DiBTagIP3D_L25MinTag4_v1");
#endif
#else
  triggerFilterList.push_back("HLT_CentralJet46_BTagIP3D_CentralJet38_BTagIP3D");
  triggerFilterList.push_back("HLT_CentralJet46_CentralJet38_DiBTagIP3D_v1");
#endif


  // chain mode
  TChain f("hbbanalysis/HBBTo4B");
  TFileCollection fc("dum","","theMergeList.txt");

  // extract generic trigger list and number of input events from ntuple files
  float lumiScaleFac = 1;
  //getHbbMetadata(fc,genericTriggerList,1000.,lumiScaleFac,_doMC);

///Igor's code
getHbbMetadata(fc,genericTriggerList,501.,lumiScaleFac,_doMC);

  f.AddFileInfoList((TCollection*) fc.GetList());
  TTree* hbbtree = &f;

  // create the btag efficiency objects

  const std::string flavlabel[nflav] = {"udsg","c","b"};
  // define ovl labels
  const std::string sovl[novl] = { "TCHPT", "TCHP6", "CSVT", "SSVHPT" };
  string btdiscr[novl] = { "TCHP", "TCHP", "CSV", "SSVHP" };
  // corresponding TCHP cuts
  double btcut[novl] = { 3.41, 6 , 0.898, 2 };
  bool bbPurityCorr[novl] = { true, true, false, false };

  // define flavor pair classes
  const int nfcDijet = 6;
  string sfcDijet[nfcDijet] = { "bb", "bc", "bq", "cc", "cq", "qq" };
  const int nfc3rd = 3;
  string sfc3rd[nfc3rd] = { "q", "c", "b" };
  // flavor triplet classes
  const int nfcTrip = 6;
  string sfcTrip[nfcTrip] = { "bbb", "bbc", "bbq", "bcb", "bqb", "non-bb" };
  // the following are used for templates
  const int nfc = 3;  
  string sfc[nfc] = { "q", "c", "b" };
  const int ncateg = 3; // position of untagged jet
  const int ncorr = 2; // correction levels
  const int ntpat = 1; // trigger pattern (online btag), set to 8 later

  // make pointers to appropriate btag discriminants
  float* theBJetTag[novl];
  for (int iovl=0; iovl<novl; ++iovl) {
    if (btdiscr[iovl] == "TCHP") {
      theBJetTag[iovl] = tcHPBJetTag;
    }  else if (btdiscr[iovl] == "CSV") {
      theBJetTag[iovl] = combSVBJetTag;
    }  else if (btdiscr[iovl] == "SSVHP") {
      theBJetTag[iovl] = svHPBJetTag;
    } else {
      std::cout << "invalid btag discriminant " << btdiscr[iovl] << std::endl;
    }
  }

  //  bTagEff* bTagEffOffline = new bTagEff("/afs/naf.desy.de/user/r/rmankel/scratch/HbbPat/CMSSW_4_2_3_patch2/src/DesyHiggsAnalyses/HBBAnalyses/test/results/job014/QCD_Pt-N14/pteta-q-nb30/plotmistag-b.root","offline",flavlabel,sovl);
  bTagEff* bTagEffOffline = new bTagEff("/afs/naf.desy.de/user/r/rmankel/scratch/HbbPat/CMSSW_4_2_3_patch2/src/DesyHiggsAnalyses/HBBAnalyses/test/results/job014/QCD_Pt-N14/btagMatrixOffline-csv/plotmistag-b.root","offline",flavlabel,sovl); 



 // this one for smp smoothing
   bTagEff* bTagReleffOnline = new bTagEff("/afs/naf.desy.de/user/r/rmankel/scratch/HbbPat/CMSSW_4_2_3_patch2/src/DesyHiggsAnalyses/HBBAnalyses/test/results/job014/QCD_Pt_bEnriched-T14b/btagMatrixOnline-csv/smooth/btagMatrixOnline-csv-smpf.root","online",flavlabel,sovl);

  // open the bbPurity correction functions
  TF1* fbbfrac[novl][ncateg];
  TFile* bbPur = new TFile("/afs/naf.desy.de/user/r/rmankel/scratch/HbbPat/CMSSW_4_2_3_patch2/src/DesyHiggsAnalyses/HBBAnalyses/test/calib/bbPurity.root");
  for (int iovl=0; iovl<novl; ++iovl) {
    for (int icateg=0; icateg<ncateg; ++icateg) {
      fbbfrac[iovl][icateg] = (TF1*) bbPur->Get( Form("fbbfrac-ww-%s-Cat%d",sovl[iovl].c_str(),icateg) );
      if ( fbbfrac[iovl][icateg] == NULL ) {
	std::cout << "bbPur correction function not found for" << sovl[iovl].c_str() << " categ " << icateg << std::endl;
	if (bbPurityCorr[iovl]) return;
      }
    }
  }

  TFile* hout = new TFile("readtribn.root","recreate");

  //const bool oldAddDir = TH1::AddDirectoryStatus();
  TH1::AddDirectory(true);

  TH1::SetDefaultSumw2();

   nentries = (Int_t) hbbtree->GetEntries();
  cout << "Number of events in ntuple: " << nentries << endl;

//#include "/afs/naf.desy.de/user/r/rmankel/scratch/HbbPat/CMSSW_4_2_3_patch2/src/DesyHiggsAnalyses/HBBAnalyses/test/readhbbset.C"
#include "readhbbset.C"



 // Histogram handler
  std::map<std::string, TH1F *> H1_;
  std::map<std::string, TH2F *> H2_;

  // counter histogram for flavor pair classes
  TH1F* hfc[novl];
  TH1F* hfcm[novl];
  TH1F* hfcww[novl];
  TH1F* hfcmww[novl];
  for (int iovl=0; iovl<novl; ++iovl) {
    hfc[iovl] = new TH1F(Form("hfc%s",sovl[iovl].c_str()),Form("Flavor triple code 3* %s",sovl[iovl].c_str()),19,-1.5,17.5);
    hfcm[iovl] = new TH1F(Form("hfcm%s",sovl[iovl].c_str()),Form("Flavor triple code 3* %s, mass 100-140",sovl[iovl].c_str()),19,-1.5,17.5);
    hfcww[iovl] = new TH1F(Form("hfc%sww",sovl[iovl].c_str()),Form("Flavor triple code weighted as 3* %s",sovl[iovl].c_str()),19,-1.5,17.5);
    hfcmww[iovl] = new TH1F(Form("hfcm%sww",sovl[iovl].c_str()),Form("Flavor triple code weighted as 3* %s, mass 100-140",sovl[iovl].c_str()),19,-1.5,17.5);
    for (int ifcDijet=0; ifcDijet<nfcDijet; ++ifcDijet) {
      for (int ifc3rd=0; ifc3rd<nfc3rd; ++ifc3rd) {
	hfc[iovl]->GetXaxis()->SetBinLabel(nfcDijet* ifc3rd + ifcDijet +2,Form("%s%s",sfcDijet[ifcDijet].c_str(),sfc3rd[ifc3rd].c_str()));
	hfcm[iovl]->GetXaxis()->SetBinLabel(nfcDijet* ifc3rd + ifcDijet +2,Form("%s%s",sfcDijet[ifcDijet].c_str(),sfc3rd[ifc3rd].c_str()));
	hfcww[iovl]->GetXaxis()->SetBinLabel(nfcDijet* ifc3rd + ifcDijet +2,Form("%s%s",sfcDijet[ifcDijet].c_str(),sfc3rd[ifc3rd].c_str()));
	hfcmww[iovl]->GetXaxis()->SetBinLabel(nfcDijet* ifc3rd + ifcDijet +2,Form("%s%s",sfcDijet[ifcDijet].c_str(),sfc3rd[ifc3rd].c_str()));
      }
    }
  }
  TH1F* hfctrip = new TH1F("hfctrip","Flavor condensed triple code",7,-1.5,6.5);
  for (int ifcTrip=0; ifcTrip<nfcTrip; ++ifcTrip) {
    hfctrip->GetXaxis()->SetBinLabel(ifcTrip+2,sfcTrip[ifcTrip].c_str());
  }

  const int minJet = 3;

  // book histogram arrays
  TrigHistArray* amdijetbt[novl];
  TrigHistArray* amdijetww[novl];
  TrigHistArray* aptjetbt[novl][minJet];
  TrigHistArray* aptjetww[novl][minJet];
  TrigHistArray* amdijetfcbt[novl][nfcTrip];
  TrigHistArray* amdijetfcww[novl][nfcTrip];
  TrigHistArray* aptjetfcbt[novl][nfcTrip][minJet];
  TrigHistArray* aptjetfcww[novl][nfcTrip][minJet];

  for (int iovl=0; iovl<novl; ++iovl) {
    amdijetbt[iovl] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
			     Form("mjjbt%s",sovl[iovl].c_str()),
			     Form("m(Jet1Jet2) 3*%s",sovl[iovl].c_str()),50,0,500 );
    for (int ifcTrip=0; ifcTrip<nfcTrip; ++ifcTrip) {
      amdijetfcbt[iovl][ifcTrip] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
						     Form("mjjbt%sfc%s",sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),
			     Form("m(Jet1Jet2) 3*%s fc %s",sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),50,0,500 );
    }
#ifdef MC
    amdijetww[iovl] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
					Form("mjjww%s",sovl[iovl].c_str()),
					Form("m(Jet1Jet2) weighted as 3*%s",sovl[iovl].c_str()),50,0,500 );
    for (int ifcTrip=0; ifcTrip<nfcTrip; ++ifcTrip) {
      amdijetfcww[iovl][ifcTrip] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
						     Form("mjjww%sfc%s",sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),
			     Form("m(Jet1Jet2) weighted as 3*%s fc %s",sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),50,0,500 );
    }
#endif
    for (int iJ=0; iJ<minJet; ++iJ) {
      aptjetbt[iovl][iJ] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
				   Form("ptj%dbt%s",iJ,sovl[iovl].c_str()),
				   Form("pt leading jet %d 3*%s",iJ,sovl[iovl].c_str()),30,0,300);
      for (int ifcTrip=0; ifcTrip<nfcTrip; ++ifcTrip) {
	aptjetfcbt[iovl][ifcTrip][iJ] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
				   Form("ptj%dbt%sfc%s",iJ,sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),
				   Form("pt leading jet %d 3*%s fc %s",iJ,sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),30,0,300);
      }
#ifdef MC
      aptjetww[iovl][iJ] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
				   Form("ptj%dww%s",iJ,sovl[iovl].c_str()),
				   Form("pt leading jet %d weighted as 3*%s",iJ,sovl[iovl].c_str()),30,0,300);
      for (int ifcTrip=0; ifcTrip<nfcTrip; ++ifcTrip) {
	aptjetfcww[iovl][ifcTrip][iJ] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
				   Form("ptj%dww%sfc%s",iJ,sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),
				   Form("pt leading jet %d weighted as 3*%s fc %s",iJ,sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),30,0,300);
      }
#endif
    }
  }

  TrigHistArray* aptJet[3];
  TrigHistArray* aetaJet[3];
  TrigHistArray* aphiJet[3];
  TrigHistArray* tcheJet[3];
  TrigHistArray* tchpJet[3];
  TrigHistArray* aIsJetWithBtag[3];

  for (int ii=0; ii<3; ++ii) {
    aptJet[ii] = new TrigHistArray(&genericTriggerList,&triggerFilterList,Form("ptj%d",ii),Form("pt of %dth leading jet",ii+1),25,0,200);
    aetaJet[ii] = new TrigHistArray(&genericTriggerList,&triggerFilterList,Form("aetaj%d",ii),Form("eta of %dth leading jet",ii+1), 100, -3.5, 3.5);
    aphiJet[ii] = new TrigHistArray(&genericTriggerList,&triggerFilterList,Form("aphij%d",ii),Form("phi of %dth leading jet",ii+1), 100, -3.15, 3.15); 
    tcheJet[ii] = new TrigHistArray(&genericTriggerList,&triggerFilterList,Form("tchej%d",ii),Form("TCHE of %dth leading jet, ptjet>50",ii+1), 60, -10, 20 );
    tchpJet[ii] = new TrigHistArray(&genericTriggerList,&triggerFilterList,Form("tchpj%d",ii),Form("TCHP of %dth leading jet, ptjet>50",ii+1), 60, -10, 20 );
    aIsJetWithBtag[ii] = new TrigHistArray(&genericTriggerList,&triggerFilterList,Form("isWithBtagj%d",ii),Form("IsJetWithBtag for %dth leading jet",ii+1),10,-3.5,6.5);
  }

  // templates
  TrigHistArray* amTemplate[nfc][novl][ncateg][ncorr][ntpat];
  for (int ifc=0; ifc<nfc; ++ifc) {
    for (int iovl=0; iovl<novl; ++iovl) {
      for (int icateg=0; icateg<ncateg; ++icateg) {
	for (int icorr=0; icorr<ncorr; ++icorr) {
	  for (int itpat=0; itpat<ntpat; ++itpat) {
	    amTemplate[ifc][iovl][icateg][icorr][itpat]
	      = new TrigHistArray(&genericTriggerList,&triggerFilterList,
				  Form("mjjTemp_%s_%s_Cat%dCorr%dTpat%d",sfc[ifc].c_str(),sovl[iovl].c_str(),icateg,icorr,itpat),
				  Form("mJet1Jet2 Template, %s %s Cat%d Corr%d Tpat%d",sfc[ifc].c_str(),sovl[iovl].c_str(),icateg,icorr,itpat),
				  50,0,500);
	  }
	}
      }
    }
  }
  
  // statistics for online btag pattern (double btag)
  TrigHistArray* atpat[novl][ncateg];
  TrigHistArray* atpatall[novl];
  for (int iovl=0; iovl<novl; ++iovl) {
      for (int icateg=0; icateg<ncateg; ++icateg) {
	atpat[iovl][icateg] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
						Form("tpat_%s_Cat%d",sovl[iovl].c_str(),icateg),
						Form("Online btag trig pattern, 2*%s, Cat%d",sovl[iovl].c_str(),icateg),9,-0.5,8.5);
      }
      atpatall[iovl] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
					 Form("tpat_%s_all",sovl[iovl].c_str()),
					 Form("Online btag trig pattern, 2*%s, all",sovl[iovl].c_str()),9,-0.5,8.5);
  }

  // statistics for online btag pattern (triple btag)
  TrigHistArray* atpattripall[novl];
  for (int iovl=0; iovl<novl; ++iovl) {
    atpattripall[iovl] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
					   Form("tpattrip_%s_all",sovl[iovl].c_str()),
					   Form("Online btag trig pattern, 3*%s, all",sovl[iovl].c_str()),9,-0.5,8.5);
  }

#ifdef MC
  // double btag purity histograms
  TH1F* hmdibbt[nfcDijet][novl][ncateg];
  TH1F* hmdibww[nfcDijet][novl][ncateg];
  for (int ifcDijet=0; ifcDijet<nfcDijet; ++ifcDijet) {
    for (int iovl=0; iovl<novl; ++iovl) {
      for (int icateg=0; icateg<ncateg; ++icateg) {
	hmdibbt[ifcDijet][iovl][icateg] = new TH1F(Form("mdib_%s_%s_Cat%dbt",sfcDijet[ifcDijet].c_str(),sovl[iovl].c_str(),icateg),
				   Form("Jet1Jet2 %s 2* %s Cat%d",sfcDijet[ifcDijet].c_str(),sovl[iovl].c_str(),icateg),
						 50,0,500);
	hmdibww[ifcDijet][iovl][icateg] = new TH1F(Form("mdib_%s_%s_Cat%dww",sfcDijet[ifcDijet].c_str(),sovl[iovl].c_str(),icateg),
				   Form("Jet1Jet2 %s weighted as 2* %s Cat%d",sfcDijet[ifcDijet].c_str(),sovl[iovl].c_str(),icateg),
						 50,0,500);
      }
    }
  }
#endif

  TrigHistArray* aMatchPattern = new TrigHistArray(&genericTriggerList,&triggerFilterList,"matchPattern","matchPattern",20,-0.5,19.5);
  
  TrigHistArray* adphiJet = new TrigHistArray(&genericTriggerList,&triggerFilterList,"adphij","delta-phi of two leading jets", 100, -3.15, 3.15); 
  
///Igor's code here

/// _MVA_ selection here: Ini of all need support
///It's done in the external driver program run.C
//gROOT->ProcessLine(".L MVACalculator.h+");
//gROOT->ProcessLine(".L MVAComputer.h+");
///Create MVAComputers


//std::string treexml = "/afs/naf.desy.de/user/m/marfin/scratch/releases/CMSSW_4_2_3/src/AlpgenAnalysis/AlpgenAnalysis2/data/tree.xml";
//InputVariablesHistograms * MVAInput=0;
//InputVariablesHistograms * MVAInputTrig=0;


std::string treexml = "/data/user/marfin/MVAPythia/tree0.xml";
MVATrainer * trainer0=0; ///before weighting 
trainer0 = new MVATrainer("MVATrainer0","MVATrainer0",treexml);

treexml = "/data/user/marfin/MVAPythia/tree1.xml";
MVATrainer * trainer1=0; /// only lumi-weight   (changed!! only first jet btag-offline weighting)
trainer1 = new MVATrainer("MVATrainer1","MVATrainer1",treexml);

treexml = "/data/user/marfin/MVAPythia/tree2.xml";
MVATrainer * trainer2=0; ///triple btag-offline weighting 
//trainer2 = new MVATrainer("MVATrainer2","MVATrainer2",treexml);


treexml = "/data/user/marfin/MVAPythia/tree3.xml";
MVATrainer * trainer3=0; ///double btag-online weighting 
//trainer3 = new MVATrainer("MVATrainer3","MVATrainer3",treexml);

treexml = "/data/user/marfin/MVAPythia/tree4.xml";
MVATrainer * trainer4=0; ///final weight 
trainer4 = new MVATrainer("MVATrainer4","MVATrainer4",treexml);

treexml = "/data/user/marfin/MVAPythia/tree44.xml";
MVATrainer * trainer44=0; ///final weight without lumi weight
//trainer44 = new MVATrainer("MVATrainer44","MVATrainer44",treexml);


//MVAInput = new InputVariablesHistograms("Input","Input",mvaPath,mvaMethods);
//MVAInputTrig = new InputVariablesHistograms("InputTrig","InputTrig",mvaPath,mvaMethods);

MVAComputer2 * comp=0;

#ifdef MVA
std::string mvaPath = "/data/user/marfin/MVAPythia/tmva/weights";
//mvaPath = "/data/user/marfin/MVAPythia/tmva/weights_ptTHCP6_democratic";
mvaPath = "/data/user/marfin/MVAPythia/tmva/weights";
std::vector<std::string>  mvaMethods;
mvaMethods.push_back("Likelihood"); /// add any other methods if you are sure that they leave at mvaPath
comp = new MVAComputer2("MVAComputer2","MVAComputer2",mvaPath,mvaMethods);
#endif


treexml = "/data/user/marfin/MVAPythia/tree5.xml";
MVATrainer * trainer5=0; ///to fill with mva response

#ifdef MVA
trainer5 = new MVATrainer("MVATrainer5","MVATrainer5",treexml);
#endif


treexml = "/data/user/marfin/MVAPythia/tree6.xml";
MVATrainer * trainer6=0; ///to fill after mva selection

#ifdef MVA
trainer6 = new MVATrainer("MVATrainer6","MVATrainer6",treexml);
#endif





///_MVA_ ini stop


//nentries=5e4;
if (nentries!=nInputEvents ){

//cout<<"nentries is inconsistent with nInputEvents, redefine lumiScaleFac"<<endl;
//cout<<"Old lumiScaleFac "<<lumiScaleFac<<endl;
//lumiScaleFac*=float(nentries)/float(nInputEvents);
//cout<<"New lumiScaleFac "<<lumiScaleFac<<endl;



}

#ifdef USERNEVT
nentries = USERNEVT;
#endif

  // loop over the tree
  for (Int_t i=0; i<nentries    ; i++) {
    hbbtree->GetEntry(i);


if (i%1000==0) {
//if (i%1==0) {
cout<<"Event # "<<i<<endl;
}

    // this is mainly for counting events before any analysis cuts
    
    double jetPtMin[3] = { 46, 38, 15};  // Roberval's cuts, but no extra cut on 3rd jet

#ifdef PTMAX300
    double jetPtMax[3] = {300, 300, 300}; // restrict to range where we have btag eff 
#else
    double jetPtMax[3] = {3500, 3500, 3500}; // restrict to range where we have btag eff 
#endif

    double maxEta = 2.6;

    std::vector<int> leadingJets;

///Igor's code
std::vector<math::XYZTLorentzVector> _vecJets ;
std::vector<double> * res;
double minPt=15e0; /// some cut on 

    double maxEta2 = 2.6;
    double minPt2=5e0 ;


    // compute the weight for this event
    float weight = 1;
#ifdef LUMIWEIGHTING
    weight = lumiScaleFac;
#endif

    bool matchOK = true;






//if (MVAInput) MVAInput->FillHistograms(_vecJets);



  std::vector<int> theFlav3;

    
	





    int nJet = 0;
    // loop over the jets
    for (int iJet=0; iJet<numberOfJets; ++iJet) {

bool _isplaced=false;


//      if (nJet >= minJet) break;
  	if (nJet<minJet)
      if ( (fabs(etaJet[iJet])<maxEta) )
      if ( (numberOfConstituentsInJet[iJet] > 1) ) 
      if ( (ptJet[iJet] > jetPtMin[nJet]) && (ptJet[iJet] < jetPtMax[nJet]) ) {
#ifdef MC
	theFlav3.push_back(jetFlavorCode(iJet));
#endif
	leadingJets.push_back(iJet);
         _vecJets.push_back(math::XYZTLorentzVector(pxJet[iJet], pyJet[iJet], pzJet[iJet], energyJet[iJet]));
	++nJet;
	_isplaced=true;
      }

	if (nJet == minJet && !_isplaced) 
	if ( (fabs(etaJet[iJet])<maxEta2) )
	if ( (numberOfConstituentsInJet[iJet] > 1) )
	if ( ptJet[iJet]  > minPt2 ) {
#ifdef MC
	theFlav3.push_back(jetFlavorCode(iJet));
#endif
	_vecJets.push_back(math::XYZTLorentzVector(pxJet[iJet], pyJet[iJet], pzJet[iJet], energyJet[iJet]));
	}
    }


///before any selection
if (i%10000==0)  {
if (trainer0) trainer0->ProcessEvent(_vecJets,theFlav3,weight,kTRUE);
} else { if (trainer0) trainer0->ProcessEvent(_vecJets,theFlav3,weight);}

    if (nJet < minJet) continue;

///End of preselection?


///Some check of pT oder:
if (ptJet[leadingJets[0]] < ptJet[leadingJets[1]] ) cout<<"Problem 1 with pt: "<<ptJet[leadingJets[0]] <<" "<<ptJet[leadingJets[1]]<<" event "<<  i<<endl;
if (_vecJets[0].Pt() < _vecJets[1].Pt() ) cout<<"Problem 2  with pt: "<< _vecJets[0].Pt() <<" "<<_vecJets[1].Pt()<<" event "<<  i<<endl;



    // check the matching flags
    int mJet = 0;
    int matchPat = 0;
    for (int iJet=0; iJet<numberOfJets; ++iJet) {
      if (! (ptJet[iJet] > 15) ) continue;
      if (! (fabs(etaJet[iJet])<maxEta) ) continue;
      if (! (numberOfConstituentsInJet[iJet] > 1) ) continue;
      if (isJetWithHltBtag[iJet]) {
	matchPat = matchPat | (1<<mJet);
      }
      ++mJet;
    }
    aMatchPattern->fill(trgAccept,matchPat);




#ifdef MATCHING
    matchOK = false;
    if (nJet>=minJet ) {
      matchOK = isJetWithHltBtag[leadingJets[0]] && isJetWithHltBtag[leadingJets[1]];
    }
#endif









///Igor's code here

bool _isOk=false;

float _deltaRj = -1;
    if ( (nJet>=minJet) && matchOK ) {

      float dphij =  phiJet[leadingJets[1]] - phiJet[leadingJets[0]];
      if (dphij>3.1415926) dphij -= 2*3.1415926;
      if (dphij<-3.1415926) dphij += 2*3.1415926;
      adphiJet->fill(trgAccept,dphij,weight);
      float detaj = etaJet[leadingJets[1]] - etaJet[leadingJets[0]];
      _deltaRj = sqrt( dphij*dphij + detaj*detaj );
    }


///Flavor determination

//    int theFlav2[minJet];
    std::vector<int> theFlav2;

#ifdef MC
if ( (nJet>=minJet) && matchOK ) {

    // determine jet flavor from MC info
    for (int iJ=0; iJ<minJet; ++iJ) {
//      theFlav2[iJ] = jetFlavorCode(leadingJets[iJ]);
      theFlav2.push_back(jetFlavorCode(leadingJets[iJ]));
    }


}
#endif


  // define ovl labels
 /// const std::string sovl[novl] = { "TCHPT", "TCHP6" };
//	int iiovl=1; /// "TCHP6" 
	int iiovl=0; /// "TCHPT" 


#ifdef DATA


//no cut-based btagging  in MC
//in data cut-based

   bool _btag=false;
if ( (tcHPBJetTag[leadingJets[0]]>btcut[iiovl]) && (tcHPBJetTag[leadingJets[1]]>btcut[iiovl]) &&  (tcHPBJetTag[leadingJets[2]]>btcut[iiovl]))
_btag=true;

/**
  if ( (tcHPBJetTag[leadingJets[1]]>btcut[0]) && (tcHPBJetTag[leadingJets[2]]>btcut[0]) ) {
   _btag = true;
  }
  if ((tcHPBJetTag[leadingJets[0]]>btcut[0]) && (tcHPBJetTag[leadingJets[2]]>btcut[0]) ) {
    _btag = true;
  }
  if ((tcHPBJetTag[leadingJets[0]]>btcut[0]) && (tcHPBJetTag[leadingJets[1]]>btcut[0]) ) {
    _btag = true;
  };

**/

if ( (nJet>=minJet) && (_deltaRj>1) && matchOK && _btag )  
if (trgAccept  & (1<<0) || trgAccept  & (1<<1) ) ///triggered data by Trigger0 or Trigger1
_isOk=true;

#endif



///use instead weights 

double wtot2=1e0;
double sumdb=0e0;
double wbt2[minJet];
double sumdb2=1e0;

#ifdef MC

if ( (nJet>=minJet) && (_deltaRj>1) && matchOK )  _isOk=true;


///OfflineBtag (triple-tag) weigting

          // now do the same without btag cut and with weighting
          for (int iJ=0; iJ<minJet; ++iJ) {
            wbt2[iJ] = bTagEffOffline->eff(theFlav2[iJ],sovl[iiovl].c_str(),ptJet[leadingJets[iJ]],etaJet[leadingJets[iJ]]);
 	   wtot2 *=  wbt2[iJ];         
		}

if (_isOk) {
if (i%10000==0)  {	
//if (trainer1) trainer1->ProcessEvent(_vecJets,theFlav2,wbt2[0],kTRUE);
if (trainer1) trainer1->ProcessEvent(_vecJets,theFlav2,weight,kTRUE);
	std::cout<<"BTagOff Weight of the first" <<PrintFlavor(theFlav2[0])<<" : "<<wbt2[0]<<std::endl;
	}
else {
//if (trainer1) trainer1->ProcessEvent(_vecJets,theFlav2,wbt2[0],kFALSE);
if (trainer1) trainer1->ProcessEvent(_vecJets,theFlav2,weight,kFALSE);
}

}


if(_isOk)
//if (trainer2) trainer2->ProcessEvent(_vecJets,theFlav2,weight*wtot2);
if (i%10000==0)  {

std::cout<<"BTagOff Weight of the 1st" <<PrintFlavor(theFlav2[0])<<" : "<<wbt2[0]<<std::endl;
std::cout<<"BTagOff Weight of the 2nd" <<PrintFlavor(theFlav2[1])<<" : "<<wbt2[1]<<std::endl;
std::cout<<"BTagOff Weight of the 3rd" <<PrintFlavor(theFlav2[2])<<" : "<<wbt2[2]<<std::endl;
std::cout<<"Triple Btag weight =" <<wtot2<<" again "<<wbt2[0]*wbt2[1]*wbt2[2]<<std::endl;

if (trainer2) trainer2->ProcessEvent(_vecJets,theFlav2,wtot2,kTRUE);
} else {
 if (trainer2) trainer2->ProcessEvent(_vecJets,theFlav2,wtot2,kFALSE);

}

///OnlineBtag (double-tag)

for (int iJ=0; iJ<minJet; ++iJ) 
for (int iiJ=0; iiJ<minJet; ++iiJ) 
{
	if (iiJ==iJ) continue;
	int kkk=minJet-iJ - iiJ;
///	if (kkk<0) continue;

	double _wgt = bTagReleffOnline->eff(theFlav2[iJ],sovl[iiovl].c_str(),ptJet[leadingJets[iJ]],etaJet[leadingJets[iJ]]);
	_wgt *=  bTagReleffOnline->eff(theFlav2[iiJ],sovl[iiovl].c_str(),ptJet[leadingJets[iiJ]],etaJet[leadingJets[iiJ]]);
	_wgt *= (1e0-bTagReleffOnline->eff(theFlav2[kkk],sovl[iiovl].c_str(),ptJet[leadingJets[kkk]],etaJet[leadingJets[kkk]]));

	sumdb+=_wgt;	
}


#ifdef  TESTIGORCODE
stringstream _str2;
_str2<<" online-btag weight  "<<sumdb/2e0;
COUT(_str2.str());
_str2.str("");
#endif

///	wtot2 *= sumdb/2e0;






        


if  (_isOk)
{

sumdb2=sumdb/2e0;

if (i%10000==0)  {

std::cout<<"BTagOnline Weight " <<" : "<<sumdb2<<std::endl;

if (trainer3) trainer3->ProcessEvent(_vecJets,theFlav2,sumdb2,kTRUE);
} else {
if (trainer3) trainer3->ProcessEvent(_vecJets,theFlav2,sumdb2,kFALSE);

}


 
} /// if Ok

#endif ///if MC




if (_isOk ) {
 _nevt_before_mva_trig++;

if (trainer4) trainer4->ProcessEvent(_vecJets,theFlav2,weight*wtot2*sumdb2);
if (trainer44) trainer44->ProcessEvent(_vecJets,theFlav2,wtot2*sumdb2);
}

///Here mva selection

if (_isOk)
if (comp) {

std::vector<double> * res;

comp->ProcessEvent(_vecJets);
res = comp->GetOutPut();

if ( res && res->size()>0) {
if (trainer5) trainer5->ProcessEvent(_vecJets,theFlav2,weight*wtot2*sumdb2,kFALSE,res,0); ///to fill mva response as well


#if defined(MVACUTLESS) && !defined(MVACUTMORE)

if (res->at(0) > MVACUTLESS) continue;

#endif


#if !defined(MVACUTLESS) && defined(MVACUTMORE)

if (res->at(0) < MVACUTLESS) continue;

#endif


if (_isOk) if (trainer6) trainer6->ProcessEvent(_vecJets,theFlav2,weight*wtot2*sumdb2,kFALSE,res,0); ///to fill after mva response


} /// if size()>0


} ///comp



if (_isOk ){
#if defined(MVACUTLESS) && !defined(MVACUTMORE)
if (i%10000==0) cout<<"I'm passed cut "<< MVACUTLESS<<endl;
#endif

#if !defined(MVACUTLESS) && defined(MVACUTMORE)
if (i%10000==0) cout<<"I'm passed cut "<< MVACUTMORE<<endl;
#endif


  _nevt_after_mva_trig++;

}




///_MVA_ training stop






    // fill jet histograms
    float deltaRj = -1;
    if ( (nJet>=minJet) && matchOK ) {

      for (int iJ=0; iJ<minJet; ++iJ) {
	// jet kinematics
	aptJet[iJ]->fill(trgAccept,ptJet[leadingJets[iJ]],weight);
	aetaJet[iJ]->fill(trgAccept,etaJet[leadingJets[iJ]],weight);
	aphiJet[iJ]->fill(trgAccept,phiJet[leadingJets[iJ]],weight);
	aIsJetWithBtag[iJ]->fill(trgAccept,isJetWithHltBtag[leadingJets[iJ]],weight);
      }
      float dphij =  phiJet[leadingJets[1]] - phiJet[leadingJets[0]];
      if (dphij>3.1415926) dphij -= 2*3.1415926;
      if (dphij<-3.1415926) dphij += 2*3.1415926;
      adphiJet->fill(trgAccept,dphij,weight);
      float detaj = etaJet[leadingJets[1]] - etaJet[leadingJets[0]];
      deltaRj = sqrt( dphij*dphij + detaj*detaj );
    }

    // determine jet flavor from MC info
    int theFlav[minJet];
    for (int iJ=0; iJ<minJet; ++iJ) {
      theFlav[iJ] = jetFlavorCode(leadingJets[iJ]);
    }

    // determine the dijet flavor code for the first two leading jets
    int theFcDijet = diJetFlavorCode( leadingJets[0], leadingJets[1] );

    // now combine with flavor of third jet
    int theFc3rd = theFlav[2];
    int theFc = nfcDijet * theFc3rd + theFcDijet;

    // now make the triplet code
    int theFcTrip = -1;
    switch (theFcDijet) {
    case 0:  // (bb)
      switch (theFc3rd) {
      case 2: // bbb
	theFcTrip = 0;
	break;
      case 1: // bbc
	theFcTrip = 1;
	break;
      case 0: // bbq
	theFcTrip = 2;
	break;
      default:
	std::cout << "Bad 3rd flavor code " << theFcDijet << " " << theFc3rd << std::endl;
      }
      break;
    case 1:  // (bc)
      switch (theFc3rd) {
      case 2: // bcb
	theFcTrip = 3;
	break;
      case 1:
      case 0:
	theFcTrip = 5; // non-bb
	break;
      default:
	std::cout << "Bad 3rd flavor code " << theFcDijet << " " << theFc3rd << std::endl;
      }
      break;
    case 2:  // (bq)
      switch (theFc3rd) {
      case 2: // bqb
	theFcTrip = 4;
	break;
      case 1:
      case 0:
	theFcTrip = 5; // non-bb
	break;
      }
      break;
    case 3:
    case 4:
    case 5:
      theFcTrip = 5; // non-bb
      break;
    default:
      std::cout << "Bad flavor codes " << theFcDijet << " " << theFc3rd << std::endl;
    }

    hfctrip->Fill(theFcTrip,weight);


    // compute invariant mass of two leading jets
    if ( (nJet>=minJet) && (deltaRj>1) && matchOK ) {
      float energyTot = energyJet[leadingJets[0]] + energyJet[leadingJets[1]];
      float pxTot = pxJet[leadingJets[0]] + pxJet[leadingJets[1]];
      float pyTot = pyJet[leadingJets[0]] + pyJet[leadingJets[1]];
      float pzTot = pzJet[leadingJets[0]] + pzJet[leadingJets[1]];
      
      float dijet_mass_sq = energyTot * energyTot - pxTot * pxTot - pyTot * pyTot - pzTot * pzTot;

      if (dijet_mass_sq >= 0) {
	float dijet_mass = sqrt( dijet_mass_sq );

	// fill the triple-btag histograms
	for (int iovl=0; iovl<novl; ++iovl) {
	  if ( (theBJetTag[iovl][leadingJets[0]]>btcut[iovl]) &&  (theBJetTag[iovl][leadingJets[1]]>btcut[iovl])
	       && (theBJetTag[iovl][leadingJets[2]]>btcut[iovl]) ) {
	    amdijetbt[iovl]->fill(trgAccept,dijet_mass,weight);
	    // triple-jet specific flavor code histograms
	    if (theFcTrip>=0) amdijetfcbt[iovl][theFcTrip]->fill(trgAccept,dijet_mass,weight);
	    // fc counting
	    hfc[iovl]->Fill( theFc,weight );
	    if ( (dijet_mass > 100) && (dijet_mass < 140 ) ) hfcm[iovl]->Fill( theFc,weight );
	    for (int iJ=0; iJ<minJet; ++iJ) {
	      aptjetbt[iovl][iJ]->fill(trgAccept,ptJet[leadingJets[iJ]],weight);
	      if (theFcTrip>=0) aptjetfcbt[iovl][theFcTrip][iJ]->fill(trgAccept,ptJet[leadingJets[iJ]],weight);
	    }
	    // tpat counting
	    int theTpat = 4 * int(isJetWithHltBtag[leadingJets[2]]) + 2 * int(isJetWithHltBtag[leadingJets[1]])
	      + int(isJetWithHltBtag[leadingJets[0]]);
	    atpattripall[iovl]->fill(trgAccept,float(theTpat),weight);
	  }

#ifdef MC
	  // now do the same without btag cut and with weighting
	  double wbt[minJet];
	  for (int iJ=0; iJ<minJet; ++iJ) {
	    wbt[iJ] = bTagEffOffline->eff(theFlav[iJ],sovl[iovl].c_str(),ptJet[leadingJets[iJ]],etaJet[leadingJets[iJ]]);
	  }
	  double wtot = weight * wbt[0] * wbt[1] * wbt[2];
	  amdijetww[iovl]->fill(trgAccept,dijet_mass,wtot);
	  // triple-jet specific flavor code histograms
	  if (theFcTrip>=0) amdijetfcww[iovl][theFcTrip]->fill(trgAccept,dijet_mass,wtot);
	  // fc counting
	  hfcww[iovl]->Fill( theFc,wtot );
	  if ( (dijet_mass > 100) && (dijet_mass < 140 ) ) hfcmww[iovl]->Fill( theFc,wtot );
	  for (int iJ=0; iJ<minJet; ++iJ) {
	      aptjetww[iovl][iJ]->fill(trgAccept,ptJet[leadingJets[iJ]],wtot);
	      if (theFcTrip>=0) aptjetfcww[iovl][theFcTrip][iJ]->fill(trgAccept,ptJet[leadingJets[iJ]],wtot);
	  }
#endif
	}

	// background templates

	for (int iovl=0; iovl<novl; ++iovl) {
	  // set category membership for this btag selection
	  bool lcateg[3] = {false, false, false};
	  if ( (theBJetTag[iovl][leadingJets[1]]>btcut[iovl]) && (theBJetTag[iovl][leadingJets[2]]>btcut[iovl]) ) {
	    lcateg[0] = true;
	  }
	  if ((theBJetTag[iovl][leadingJets[0]]>btcut[iovl]) && (theBJetTag[iovl][leadingJets[2]]>btcut[iovl]) ) {
	    lcateg[1] = true;
	  }
	  if ((theBJetTag[iovl][leadingJets[0]]>btcut[iovl]) && (theBJetTag[iovl][leadingJets[1]]>btcut[iovl]) ) {
	    lcateg[2] = true;
	  };
	  // determine the online btag pattern
	  int theTpat = 4 * int(isJetWithHltBtag[leadingJets[2]]) + 2 * int(isJetWithHltBtag[leadingJets[1]])
	    + int(isJetWithHltBtag[leadingJets[0]]);
	  for (int icateg=0; icateg<ncateg; ++icateg) {
	    if (lcateg[icateg]) atpat[iovl][icateg]->fill(trgAccept,float(theTpat),weight);
	  }
	  atpatall[iovl]->fill(trgAccept,float(theTpat),weight);

	  // fill the templates
	  for (int ifc=0; ifc<nfc; ++ifc) {
	    for (int icateg=0; icateg<ncateg; ++icateg) {
	      if (lcateg[icateg]) {
		float mistagWeight = bTagEffOffline->eff(ifc,sovl[iovl].c_str(),ptJet[leadingJets[icateg]],etaJet[leadingJets[icateg]]);
		for (int icorr=0; icorr<ncorr; ++icorr) {
		  double wbbpur = 1;
		  if (bbPurityCorr[iovl] && (icorr == 1)) wbbpur =  fbbfrac[iovl][icateg]->Eval( dijet_mass );
		  for (int itpat=0; itpat<ntpat; ++itpat) {
		    // no additional corrections or trigger patterns yet
		    amTemplate[ifc][iovl][icateg][icorr][itpat]->fill(trgAccept,dijet_mass,weight * mistagWeight * wbbpur);
		  }
		}
	      }
	    }
	  }
	}

#ifdef MC
	// here we fill the double btag purity histograms
	// auxiliary array to perform permutations among first three jets
	int shift[6] = { 0, 1, 2, 0, 1, 2 };
	
	for (int iovl=0; iovl<novl; ++iovl) {
	  for (int icateg=0; icateg<ncateg; ++icateg) {
	    int jetA = leadingJets[shift[icateg+1]];
	    int jetB = leadingJets[shift[icateg+2]];
	    int dibFcDijet = diJetFlavorCode( jetA, jetB );
	    if ( (theBJetTag[iovl][jetA]>btcut[iovl]) &&  (theBJetTag[iovl][jetB]>btcut[iovl]) ) {
	      hmdibbt[dibFcDijet][iovl][icateg]->Fill( dijet_mass, weight );
	    }
	    // same with btag weighting
	    hmdibww[dibFcDijet][iovl][icateg]->Fill( dijet_mass, weight
						     * bTagEffOffline->eff(theFlav[shift[icateg+1]],sovl[iovl].c_str(),ptJet[jetA],etaJet[jetA])
						     * bTagEffOffline->eff(theFlav[shift[icateg+2]],sovl[iovl].c_str(),ptJet[jetB],etaJet[jetB]) );
	  }
	}

#endif

      } else {
	std::cout << "Unphysical dijet_mass_sq = " << dijet_mass_sq << std::endl;
	//dijet_mass = -1;
      } // if (dijet_mass...)
    } // if (njet...)
  }


///Igor's code

stringstream _str;

_str<<" MVA Total # of events  "<<nentries;
COUT(_str.str());
_str.str("");

_str<<" MVA  #  Events come to MVA Train   "<<_nevt_before_mva_trig;
COUT(_str.str());
_str.str("");

_str<<"  # MVA # Triggered Events come after MVA    "<<_nevt_after_mva_trig;
COUT(_str.str());
_str.str("");






  // termination
  hout->Write();
  hout->Close();

if (trainer0) delete trainer0;
if (trainer1) delete trainer1;
if (trainer2) delete trainer2;
if (trainer3) delete trainer3;
if (trainer4) delete trainer4;
if (trainer44) delete trainer44;
if (trainer5) delete trainer5;
if (trainer6) delete trainer6;
if (comp) delete comp;

/*
if (MVAInputTrig) {
	delete MVAInputTrig;
}
if (MVAInput) {
delete MVAInput;
}
*/

}
